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ABSTRACT 


Two  method*  have  been  outlined  in  detail,  and  one  of  them 
has  been  mechanized,  for  calculating  acoustic  ray  paths  emanating 
from  any  point  in  a  non-uniform  transonic  flow  field  surrounding 
a  wing.  It  gives  the  ray  path,  and  the  time,  for  the  minimum  time 
of  travel  from  the  acoustic  source  point  to  the  field  point.  The 
resulting  velocity  potential  is  also  computed. 

It  was  necessary  to  establish  an  accurate  representation  of 
the  flov  characteristics  in  the  field  surrounding  the  wing.  Some 
ray  lines  travel  over  the  planform  and  into  the  surrounding  flov 
field.  It  vas  established  that  once  off  the  planform  they  do  not 
return. 

Available  methods  predict  phase  lags  based  on  the  assumption 
that  acoustic  rays  travel  in  straight  lines.  The  results  of  this 
study  show  this  to  be  a  very  poor  approximation  at  transonic  speeds. 
Therefore,  it  is  recoasended  that  the  method  presented  in  this 
report  be  fully  developed  for  the  purpose  of  calculating  generalized 
forces  on  wings  in  harmonic  motion  at  transonic  speeds.  A  computer 
program  that  would  predict  these  phase  lags  with  reasonable 
accuracy,  and  the  corresponding  flutter  characteristics  and  unsteady 
aerodynamic  loads  on  a  wing  responding  to  externally  applied  forces, 
such  as  gusts,  would  fill  an  important  gap  in  the  available 
technology . 
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INTRODUCTION 


When  an  airfoil  travels  through  the  air  at  speeds  near  the  speed 
of  sound,  the  local  speed  of  flow  varies  from  subsonic  near  the  for¬ 
ward  edges  to  supersonic  near  the  trailing  edges.  These  wide  variations 
of  speed  from  that  of  the  free- stream  characterize  the  non-uniform 
transonic  flow.  This  non- uniformity  of  the  flow  field  must  be  accounted 
for  in  accurate  calculations  of  unsteady  pressures  and  forces;  partic¬ 
ularly  their  phase  lags. 

In  order  to  determine  an  unsteady  transonic  flow  field  one  requires 
solutions  for  singularities  immersed  in  a  non-uniform  steady  flow, 
(Reference  l).  Source  solutions  for  a  mean  flow  that  varied  in  the 
x-direction  only  were  given  in  the  high-frequency  limit  by  Landahl 
(Reference  2).  Rodemich  (Reference  3)  presented  a  "box”  solution, 
based  on  pulsating  doublets,  which  assumes  a  uniform  mean  flow  at  Mach 
number  1.0.  No  exact  solutions  for  the  case  of  a  mean  flow  with 
arbitrary  spatial  variations  have  been  found,  thus  far,  but  Landahl 
proposed  the  basic  form  of  a  solution  which  removes  most  of  the  limita¬ 
tions  and  restrictions  of  these  approximate  solutions.  The  method 
focuses  attention  on  the  time  of  transmission  of  an  acoustic  signal 
from  a  pulsating  sending  source  to  a  distant  receiving  point.  The 
signal  travels  through  a  nearly  sonic  flow  field  where  the  Mach  number 
varies  in  a  prescribed  manner. 

This  report  contains  a  difference  equation  approach,  and 
differential  equation  approach  to  confuting  the  paths  and  the  trans¬ 
mission  times  for  acoustic  signals.  'Die  Independent  variable  in  the 
latter  approach  is  a  spatial  rather  than  a  time  variable.  A  pro¬ 
cedure  that  could  be  used  to  calculate  the  velocity  potentials  and 
generalized  forces  on  an  oscillating  surface  is  described. 
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POTENTIAL  OF  A  UNIT  SOURCE 


The  basic  expressions  proposed  by  Landahl  for  the  velocity  potential 
at  the  point  (x,y,z)  due  to  a  pulsating  source  at  (x0>yQz0)  are: 

(a)  for  a  source  In  a  locally  subsonic  flow  region 

0  ■  exp  {io>Tt-g(x,y,z,x  ,y  ,z  )])  (l) 

4ttR  o  o  o 

where  _ 

R  -  7(x-xo)2  +  [l-^(x,y,z)][(y-yQ)2  +  (z-zq)2] 

M  *  Local  Mach  Number 

x0>yQ>z0  “  Location  of  source  point 

g(x,y,z,x  ,y  ,z  )  ■  Time  required  for  a  disturbance  to  travel  from 

(x0,y0,*0)  t0  (*>*>*)• 

(b)  for  a  source  in  a  locally  supersonic  flow  region 

0  -  “  {exp[io>(t-g  )]  +  exp[ixn(t-g  )]1  (2) 

UnR  a  r 

where 

6  *6  (x,y,z,x  ,y  ,z  )  ■  Time  required  for  the  advancing,  reced- 

a,r  a,r  o  o  o  wave  trave^  frOB  (x  fy  >z  )  to 

(\  o  c  o 

x,y,z) 

It  is  likely  that  good  accuracy  may  be  obtained  with  use  of  the  value 
of  ga  for  uniform  flow  (in  the  siq>er sonic  case,  and  also  for  the  advancing 
wave  portion  in  the  subsonic  case).  However,  our  purpose  is  to  produce  a 
general  solution  for  g  which  applies  to  both  the  advancing  and  the  receding 
portions  of  the  wave  and  compare  values  with  those  for  uniform  flow. 

Since  the  primary  interest  is  in  wing  flows,  we  consider  that  both  the 
source  and  receiver  points  lie  in  the  x,  y-plane,  so  that  z  »  zQ  ■  0. 
Furthermore,  we  consider  that  signals  do  not  return  to  the  plane  once  they 
leave.  The  problem  is  thus  simplified  to  one  in  two  spatial  dimensions. 

Its  solution  should  be  applicable  to  a  wide  variety  of  nearly  planar  lift¬ 
ing  surfaces. 

Consider  a  signal  emanating  from  a  source  at  the  point  (xo,y0)  on  a 
wing.  A  second  point  past  which  the  signal  travels  is  located  an  incre¬ 
mental  distance  (dx,  dy)  away.  There  are  two  components  of  velocity  of 
the  signal,  a  radial  component,  C,  where  C  is  the  local  speed  of  sound  and 
an  x-component,  U,  where  U  is  the  local  speed  of  flow  over  the  wing.  A  is 
the  angle  the  radial  component  makes  with  the  negative  extension  of  the 
x-ax.  s.  The  path  of  this  wavefront  point  will  be  referred  to  as  a  "ray". 

The  shape  of  any  ray  depends  on  the  iritial  choice  of  A;  for  a  given  A,  dx 
and  dy  are  components  of  the  first  element  of  this  particular  ray  emanating 
from  (xo,y0).  The  situation  depicted  is  general  in  that  it  applies  not  only 
at  the  source,  but  at  any  point  on  the  ray  path.  Thus,  the 
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velocity  at  any  point  on  the  path  is  a  function  of  three  spatial  para¬ 
meters  which  vary  with  position,  U,  C,  and  A.  From  the  sketch,  it  is 
clear  that 

dx  =  [U(x,  y)  -  c(x,  y)  cos  A]  dt  (3) 

dy  *  C(x,  y)  dt  sin  A 

Equations  were  developed  for  two  methods  of  tracing  the  ray  path  to 
establish  the  magnitude  and  the  phase  relationship  at  field  points  to  a 
unit  source.  These  methods  are:  (l)  a  difference  equation  method,  and 
(2)  a  non-linear  differential  equation  method. 


Difference  Equation  Method 

In  this  method,  time  is  the  independent  variable.  Equations  (3)  are 
two  of  the  three  equations  needed  to  establish  the  variation  of  x,  y,  and 
A  with  time.  The  third  equation  is  obtained  by  considering  the  accelera¬ 
tion  of  the  ray  in  the  non-uniform  flow  field  (see  Figure  l). 
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Figure  1.  Velocity  Components  of  a  Sonic  Ray  Line 
In  A  Moving  Airstream 


In  term?  of  conQ>onent8  in  the  directions  of  the  rotating  unit  vectors  ^ 
and  f 

Pt,  ‘  (Crsi*-4)f'*(C-V  l)l'  (l( 

and  K,  =  {irsmsi.  *CU)a4-(C -O  t 

It  is  necessary  to  express  the  angular  velocity  A  in  terms  of  space  vari¬ 
ables.  To  do  this,  t  cor  «**  '•cr  that  at  time  t  a  second  ray  point  is  located 
at  SR  *  ,  where  BR  is  small,  and  it's  direction  of  travel  is 

/jft  *  /R,  t  .  Let  the  superscripts  (o)  and  (l)  denote  times  t  and 

tn  (■  t  +  At) .  Then  at  time  t,  0 

*  0  1 
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and 


<'J  =  /*/"  t  *:•>*< 

Subtracting  the  first  equation  from  the  second 

5/R  0>  *  (5) 

where  S0i  *  0fx 

Recalling  that  the  cross  product  of  two  vectors  Is  a  vector  normal  to  the 
plane  defined  oy  the  two  vectors,  and  has  a  Magnitude  equal  to  the  product 
of  the  two  ir  gnltudes  tines  the  sine  of  the  angle  between  them,  then 

/•  V  *  IJR0)  x  4 '(-  rxft>  rx'°  5/A  *  a)  (6) 

which  has  the  correct  sense.  When  AA  is  snail,  and  when  Equation  (5)  Is 
substituted  into  the  left  side  of  Equation  (6),  we  get 

r  t' aA) 

This  may  be  rewritten  as 

AA  _  l«*-u  «0»a) 

A*  ~  ~  ffi*" 
and  in  the  Halt  as  At  -•  0 

JL  *  -  (7) 

where  the  operator  is 

A'W  *  (  a  />*  fa  +  cos  A 

and  operates  only  on  C  and  U. 

Equation  (7)  has  a  revealing  physical  interpretation.  From  Figure  4  we 
see  that  the  gradient  of  the  speed  of  sound  C,  on  forward  portions  of  the 
wing,  is  a  vector  pointing  forward  and  slightly  outward  from  the  center- 
line;  whereas,  from  Figure  3  we  see  that  the  gradient  of  the  local  flow 
speed  U  is  nearly  in  the  opposite  direction.  Although  it  is  not  apparent 
from  the  figures  because  they  are  plotted  to  different  scales,  the  magnitude 
of  the  gradient  of  U  is  about  five  times  that  of  the  gradient  of  C.  From 

the  energy  equation  <?  ♦  &  *  constant,  vU  -  -5.0  7C.  The  local  Mach 

number  is  increasing  in  the  downstream  direction.  Figure  2  shows  that, 
under  these  conditions  there  are  only  two  stable  ray  angles;  those  for 
which  the  gradient  of  C  -  U  cos  A  is  zero.  As  the  ray  propagates  through 
the  flow  field  it  will  always  tend  towards  one  of  these  two  orientations. 
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Figure  2.  Stability  of  Ray  Angles  When  The  Gradient  of 
Local  Flow  Speed  Exceeds  the  Gradient  of 
Local  Speed  of  Sound. 


We  now  write  Equations  (3)  and  (7)  in  difference  fora 


Ax  * 

b 

-  C  cos  A]  Ag 

(8-a) 

Ay  - 

b 

sin  AjAg 

(8-b) 

AA  - 

■l 

8lnA 

(8-c) 

*e  oj  ** 

where  represents  an  increment  in  disturbance  travel  time  g,  defined 
previously.  To  determine  cp(x,  y,  0,  Yo>  °)  It  1®  necessary  to  know 
a  steady  state  distribution  of  C(x,  y),  U(x,  y),  and  their  derivatives 
at  any  point  in  the  flow  field  over  the  wing  and  in  the  surrounding  flow 
field  in  the  plane  of  the  wing.  A  means  for  establishing  these  is  given 
in  Section  5«  Assume  they  are  known.  Then  the  procedure  used  is  as 
follows : 

1.  Select  any  source  point,  on  or  off  the  wing,  (xQ,  yQ). 
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2. 


Select  «  eerie*  of  initial  ray  angles,  A^ ,  1  ■  1,  2, - . 


3. 

4. 


8elect  an  Initial  Inc 


>nt  In  disturbance  travel  time,  Ag  . 

o 


For  each  of  the  raj  angles  store  x^,  y^,  sin  cos 


and  Ag 

a 


(i) 


At  x^,  y^  compute  and  store  x^  ■  x^  ♦  Ax^/2  and 
y(i)  +  y(l)  +  holding  A  constant. 

b.  Iterate  on  .  ,<‘>  .  *<‘>/2(  yM  .  7d)  ,  „d)/2> 

and  AA(x^^\  y2 ' 1  ^ )  until  they  converge  or  exceed  ten  trials. 

In  the  latter  case  replace  Ag^  bj  Ag  *V2  ud  repeat  the 
Iteration.  If  they  converge  In  three  trials  or  less,  replace 

Ag*1*  by  2Ag^. 

c.  Replace  x^  by  Xg^,  F^  >  •nd  r*turn  *>  »• 

the  solutions  presented  above  are  believed  to  be  good  approximations 
to  the  exact  solutions  for  the  following  reasons: 

1.  For  the  case  of  a  uniform  flow  they  reduce  to  the  proper  linearised 
expressions. 


2.  The  phase  of  the  disturbance  will  be  exact,  although  the  amplitude 
nay  be  slightly  In  error. 


3.  In  an  Inner  region  In  the  Immediate  neighborhood  of  the  source 

Location  (x  ,  y  ,  •  )  they  approach  the  correct  solution. 

0  0  0 

4.  For  a  one- dimensional  mean  flow  with  Mj  approaching  unity  they  re¬ 
duce  to  Landahl's  earlier  solution  (Reference  2). 

5.  In  the  limit  of  steady  flow  (®  -  0),  the  solutions  give  results 
equivalent  to  the  local  linearisation  method  of  Sprelter  and  Alksne 
(Reference  4).  Ails  has  been  demonstrated  by  Rubbert  (Reference  5). 

6.  Inamroch  as  the  proposed  approximation  only  affects  the  receding 
part  of  the  solution,  the  proper  limiting  solution  for  high  fre¬ 
quencies  (Reference  l),  should  always  be  obtained  since  then  reced¬ 
ing-wav#  effects  are  largely  cancelled  out  due  to  the  rapid  phase 
variations. 


This  method  gives  reasonable  results,  l.e.,  reasonable  based  on  a 
comparison  with  results  obtained  from  the  differential  Equation  method. 
However,  the  ray  paths  did  not  conclusively  show  the  existence  of  the 
focal  point  that  the  second  o«thod  revealed. 
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Non-Linear  Differential  Equation  Method 

From  Equations  (3)  we  may  write  the  slope  of  the  ray  path 

dd  _  Af- 
dy  sin  -A. 

and  solving  this  equation  for  cos  A,  we  get 


cos  A 


Af  ±  r 


t+  r 


(9) 


(10) 


where  _  <t£. 


The  transmission  time  from  source  to  receiving  point  is  given  by 

<u> 


where  the  Integration  is  taken  along  the  path  and 
ds  ■  /TT  r^  1  dy 

The  velocity  along  the  path  is  obtained  from  the  vector  6um  of  the  two 
velocity  con?>onents 


(12) 


have: 


V.c  (7  +  1  -  2M  cos  A  (13) 

Substituting  equations  (12),  (13) ,  and  (10)  into  equation  (ll)  we 


r 


/hTr^  d* 


/  -  zAf[  -  ] 


which  reduces  to 


( /+  rl)  d * 


*  I  C  /Aftrt^f  ZAfr/r'+t^A/^  +  r  x+/-a*  x ' 

The  radioand  in  the  denominator  is  a  perfect  square.  Thus, 


(14) 


(J  +  r*)dj* L 


C  [Afr  ?  /r  x+i-n  *] 


which  reduces  to 


7ir±  /r2*/-#1-' 

C  (#*-/) 


d> 


(15) 
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At  this  point  we  relate  the  local  acoustic  velocity,  C  •  c(x,  y), 
to  the  local  Mach  number  by  imposing  the  condition  of  conservation  of 
energy.  For  non-viscous  flow,  the  total  tenqjerature  io  conserved.  It 
is  easily  verified,  that  under  this  condition 


S' 

S~+*J 


(16) 


where  I4  -  1.4,  for  a  diatomic  gas,  has  been  used.  Substituting  Equation 
(l6)  into  Equation  (15)>  we  get 


Ys+m  r  i  /r  *’] 


(**-') 


where  the  ypper  sifpi  applies  to  receding  waves  and  the  lower  sign  to  ad¬ 
vancing  waves.  Equation  (17)  contains  all  the  elements  for  the  solution. 
However,  the  integrand  is  a  function  of  x,  y,  and  dx/dy.  This  equation 
may  be  written  in  symbolic  form 

/"*'  ,  ... 


which  suggests  the  use  of  Euler's  equation  to  find  the  minimum  time  g,  for 
the  disturbance  to  travel  to  a  field  point  (xi»yi) 

A-  &£■  -  *  o  /l8% 

df  dr  d*  (10) 

In  order  to  simplify  the  notation,  we  set 
J  (Sfrta) 
r  "  m  t-/  _ 

where  6  -  6(x,  y)  ■  / ^  +  W 


0  (ac,  y,  r) 


and  r  has  been  previously  defined.  We  will  need 


=  _ _ 

d'K  * 
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Ihen,  making  use  of  the  relationships 


dAi 

d* 


-  r 

%  7 


%  -■  flr% 


iL 

d9 


-t[ 


r#st^  a 


sol  ring  for  dr/dy,  and  combining  terms,  we  get 

ir,  _  /  fr-#(*+'ir  fi(7»l+snri 

‘  #•-/  f  **•'  J  r 


y-jlz M(A l+r)r±p ( 7* **>*jjJ*f  (rt+<')J'f/t 


(19) 


Equation  (19)  is  a  second  order,  second  degree  differential  equation 
of  the  form 

*«*,*.%) 

It  is  second  degree  because  8  represents  a  radical.  However,  it  can 
be  solved  numerically  by  any  of  the  standard  repetitive  processes.  We 
enployed  a  fourth  order  Runge-Kutta  procedure. 


There  are  certain  difficulties  that  arise  in  the  numerical  valua¬ 
tion  of  Equation  (19).  These  are  first  listed  and  interpreted  a  then 
equations  used  to  surmount  them  are  presented. 


(1)  Along  some  ray  paths  dx/dy  becomes  infinite  even  when  the  Mach 
number  is  not  equal  to  one. 

(2)  Equation  (19)  is  singular  at  Mach  number  =  1.0. 

(3)  In  the  supersonic  region,  signals  sometimes  become  trapped  on  the 
local  Mach  line.  This  happens  wh»n  cos  A  =  l/M.  Signals  tend  to 
gravitate  to  this  condition.  Such  trapped  signals  cannot  then  cross 
the  sonic  line.  They  approach  the  sonic  line  as  a  limit,  and  are 
cancelled  out  there. 


To  overcome  the  difficulty  listed  in  Item  (l),  it  is  necessary  to  U6e 
x  instead  of  y  as  the  independent  variable.  This  is  done  by  applying  the 
equation 


-/ 

(%r d*1 


(20) 
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It  is  convenient  here  to  Introduce  erne  new  notation.  Re*  write 
equation  (19)  in  the  fora 


-r^r  l'*+ *#(*%*)*.'* ( 7*W) 4- \*L 

4B(B  J  ?  (21) 

*  -*-{  *'*+*]"• 

where  the  new  notation,  together  with  sobm  other  notation  which  will  be 
ueed  later,  la  defined  aa  followa 


*'•-**■ 

/4  r 

„  --^2 

y'*& 

/?,  *  fa'-C#*-')1 

5  •  /¥ -/ 

~  =  ^2 

*1  *  //-  y'V^V; 

(22) 


(23-*) 


(23-b) 


The  Halting  fora  of  Equation  (20)  at  M  ■  1  la: 

In  the  eiqteraonlc  region,  when  the  signal  la  trapped  on  the  local  Mach 
line,  and 

co*-st *  if  ,  1%'hf 

equation  (20)  reduce  a  to 

*>*  • 


Substituting  Equation  (20)  Into  Equation  (21),  ve  get 

#  "*  ^  W-  *#(*  *  'V  *  * 

-*£-(<> 9 1+')  ;  i*0 
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A  complete  set  of  equations,  together  with  their  areas  of  applica¬ 
bility,  will  now  be  outlined. 

Corqplete  Set  of  Equations  where  Y  is  the  Independent  Variable 


- 

/ 

AB 

+ 

dt 

i 

t r (Al  ±  #,) 

d *  ‘ 

a 

At1-/ 

— 

*'V*V  fy 

> 


«L. 


'M-he 


*"l  =  + 

'/* 7«f»  '  *  7 


At  fc+Atl_ 


(2k) 

(25) 

(26) 

(27) 

(28) 

129) 


A  conplete  set  of  equations  were  alBo  developed  using  x  as  the  indepen¬ 
dent  variable.  However,  for  the  sake  of  brevity,  and  since  they  are 
obtained  by  a  single  change  of  variable,  they  will  not  be  listed  here. 
Equations  (26)  and  (27)  apply  where  an  advancing  ray  path  crosses  the 
sonic  line,  and  equations  (28),  (29)  apply  where  a  ray  path,  in  the  super¬ 
sonic  region,  becomes  trapped  on  the  local  Mach  line.  It  remains  to 
describe  the  regions  of  applicability  of  the  upper  and  lower  signs  of 
equations  (24)  and  (25).  In  what  follows,  "right  branch"  will  be  speci¬ 
fied  where  rf)  and  left  branch  will  be  specified  if  (-rrt  *  o  ). 

Here  is  the  local  value  along  the  ray  path.  The  end  points  are  not 
specified  because  for  these  points  we  use  x  as  the  independent  variable. 

The  upper  sign  is  used  for 

(1)  Subsonic,  left  branch 

(2)  Svqpersonic,  receding,  right  branch 

(3)  Supersonic,  advancing,  left  branch 

11 


Ifce  lover  sign  is  used  for 

(1)  Subsonic,  right  branch 

(2)  Supersonic,  receding,  left  branch 

(3)  Supersonic ,  advancing,  right  branch 
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TEE  NON -UN If' OHM  FLOW  FIELD 


In  the  application  of  each  of  the  methods  contained  in  this 
report,  it  is  necessary  to  know  certain  of  the  properties  of  the 
transonic  flow  field  on,  and  in  the  neighborhood  of,  the  win^. 
Figures  3  and  4  show  the  distributions  of  local  flow  speeds  and 
sonic  speeds  over  a  65°  delta,  wing  model  in  a  wind  tunnel  in  wtii-'h 
the  Mach  number  was  1.04  (taken  from  Reference  6).  Speeds  were 
computed  from  steady  state  pressure  data  at  2J  points  on  the  win,;. 
The  figures  are  intended  only  to  show  the  general  characteristics 
of  the  flow,  such  as:  (l)  The  local  sonic  line  shlftB  aft  with 
distance  from  the  centerline  but  crosses  the  leading  edge  inboard 
of  the  tip,  (2)  Mach  number  variations  in  both  the  streamwise  and 
spanvi6e  directions  must  be  considered  and  cannot  be  considered  to 
be  linear,  and  (3)  Separated  flow  is  indicated  over  the  aft  and 
Inboard  portion  of  the  wing.  To  consider  the  last  of  these  charac' 
teristics  is  beyond  the  scope  of  this  study.  However,  the  first 
two  are  amenable  to  analysis  using  available  theories  and 
techniques. 
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Figure  ».  Sonic  Speed  Distribution  on  a  65"  A  at  a  Transonic  Speed 


Mach  number  distributions  over  areas  off  the  wing  were  computed 
from  an  approximate  theoretical  solution  of  the  flow  field  that  matched 
pressure  distributions  on  the  wing.  In  order  to  avoid  a  discontinuity 
at  the  Juncture  of  the  two  regions,  a  small  transition  region  was 
defineu  over  which  the  two  functions  were  Joined  by  a  numerical  smooth¬ 
ing  technique. 

Let: 

■  Mj  (x,y)  °  Mach  number 

(x,y)  ■  Perturbation  notential 

c*  _ 

I  S  T^c(x,y)  -  Sickness  ratio 

For  a  steady- state,  non-lifting  flow 

4  *»*  =  ° 


Where  0  a  function  describing  the  variation  of  the  surface 

from  the  swan. 

Using  parametric  differentiation  with  respect  to'f,  (Reference  S) , 

$=  JCM)  -  §£ 

Equation  (30)  becomes: 

£ji>  *  hi  =  ° 

1>e*)  =  *i*C*.l) 


After  having  obtained  the  solution  of  equation  (32),  the  local  Mach  number 
distribution  is  obtained  by  relating  local  Mach  number  to  the  coefficient 
of  pressure,  (Cp).  Starting  with  the  following  basic  relations: 


Let 

r i 

_ 

Vi  -  lb. 

then 

"K. 

— 

(33) 

(t 

+ 

fz=  Ct-nsiti-nL 

(3U) 

where 

q 

ms 

Hm  &t  infinity 

q 

■i 

£/■*■><-)  elsewhere 

a 

m 

speed  of  sound 
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We  have: 


oi.  +±(f-i)v£  *  a.l+i.C.r-OviC'**-)'’ 
t-l  =  *i„  -  o' )  **?  *- 

using  equation  (  ^ 

a'  =  L1 

T’he  coefficient  of  pressure,  is  of  onler  f.l^.  and  M  is  Ofl.). 
Therefore,  to  sufficient  accuracy. 

a.L  =  a.„[i  +  ±(r-i)M£c/> J 

ZC  -  Vo.  t&>  ('-?<» 

and  from  these  relations: 

M  =  Moo  (I  -  ecr) 

!  +  C  0  p 

Noting  again  the  order  of  and  Op,  to  sufficient  a^curacv. 

or  Mt  S  /*/«,  [;  -  -£±i-  <>  ] 

’•’ouation  ( ?5)  is  the  expression  that  was  used  to  relate  local  Mach 
number  to  on  regions  off  the  vine. 

A  solution  of  eouation  .  using  the  results  of  equation  ( was 
worked  out  for  a  snecial  configuration.  The  snecial  wine  configuration  is 
depleted  In  figure  (s). 

u.  j 


Topical  Section 


Fie.  5.  A  Thin  Wing  In  Rectilineal  Flight 
The  solution  Is: 

C P(*. 3) -Crfa, °)  -  -2f[l#'S|  +l)*sl-ilsl  J 


(161 


where  f/(z)  ia  a  ster  function 
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(37) 


[-  If  S’[V  /,  -  5,/  ■  -  J /</-*>-*) 

[' b's>l<+  k) 


%&;=-jfe[l9-slt  '+■  Is+Sl£  J 

+iflfl[ia-s.P+lg**rfr'JH<z-‘t) 


X  r  — 

J  3 /VA- 


CCS1*- 


>  t  = 


J2.7T  dL  fio*  /L 


(18) 


X  -  ,  £,  -  -r^TTZ^C  1  s,  z 

x  r  £2Lii  >  * «. =  atr^  >  sz 

J  ».  J/V  A  »_  *• 


After  determining  a  distribution  of  C^,  and  its  derivatives  from  equations 
(16),  ( 37) ,  and  ( 38) .  the  Mach  number  distribution,  with  its  derivatives,  is 
computed  from  equation  (15). 
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DESCRIPTION  OF  THE  COMPUTER  PROGRAM 


The  equations  for  the  ray  paths  are  solved  in  the  following  manner: 
Let  the  independent  variable  be  y  and 


Vl 

vj 

VJ 


% 

t 


Then 


£*•  v, 

7?  -  UV..V,,?) 


These  three  simultaneous  differential  equations  are  solved  in  a  step-by- 
step  manner  by  use  of  a  standard  "SHARE"  subroutine  which  is  based  on  the 
Runge  Kutta  method.  When  dx/dy  becomes  greater  than  one,  a  variable 
change  takes  place  in  the  program,  and  x  becomes  the  independent  variable. 


A  signal  (in  the  supersonic  region)  i6  considered  "trapped"  on  the 
local  Mach  line  when 

When,  for  thia  trapped  signal,  (M-l)  <  E2,  the  integration  stops  and  a  new 
ray  line  is  started.  This  logical  flow  is  shown  in  the  chart  on  page  21. 


The  values  of^tm  used  in  the  program  are  dexermined  by  the  parameter 
(NLA).  If  (NLA)  is  an  odd  integer,  it  will  be  rounded  down  in  the  program 
to  an  even  integer.  Values  of  vary  frcn  zero  to  7f  and  from  zero  to  -/T 
in  an  arithmetic  progression. 

Computation  of  a  ray  path  (other  than  for  a  "trapped  signal")  ceases 
under  the  following  conditions: 

*  *  O 
/  ** 

|  YMAX  |  6  If  I 
A/ At  A  A  A/CA/r 


where  NCNT  is  the  number  of  points  on  the  ray  path  already  computed.  This 
logical  flow  is  shown  in  the  chart  on  page  22. 

Subroutine  DERIV  computes  the  appropriate  derivatives. 

Subroutine  CNTRL  accomplishes  variable  changes,  stores  local  values 
in  appropriate  locations  for  later  printing,  and  performs  exit  tests. 

Subroutine  FMACH  computes  the  local  Mach  number  and  the  partial 
derivatives  of  the  Mac'  number. 
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Subroutine  SONK  compute*  coordinate*  on  the  planform  where  M  -  1. 

Sample  data  sheets  with  numbers  which  have  been  used  in  a  computer 
run  are  in  Appendix  II.  The  output  sheets  are  included.  The  output 
format  is  self- exp lar  ‘‘ory,  with  the  exceptions  of  certain  test  words 
that  are  printed  ou  t  the  beginning  of  the  plots  for  each  ray- path. 

Definitions  for  the  words  can  be  found  in  the  consent  statements  at 

the  beginning  of  the  listing  in  Appendix  I.  The  values  listed  for  these 
test  words  apply  to  the  last  point  plotted  for  the  ray-path. 
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MAIN  PROGRAM 


Subroutine  SONK  Computes  Sonic  Line 


Subroutine  LIMITI  Sets  Plotting  Grid 
Limits 


Subroutine  GRAPH  Produces  Cathode 
Ray  Tube  Plots 


Subroutine  FMACH  Computes  Mach  No. 

Determines  Whether  X  or  Y  is 
Independent 

Determines  Left  or  Right  Branch 
Determines  Type  of  Source 


Runge  Kutta  Integrating  Subroutine 

Subroutine  POT  Computes  Velocity 
Potential  Along  Path  due  to  Source  at 
(Xo,  Yo) 
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SUBROUTINE  RKS3 


SUBROUTINE  RKZNT 


STORAGE 
I  ALLOCATION 


I 


INITIAL  VALUES 

OF  VARIABLES  — 1 


CALL  DERIV  COMPUTE 
INITIAL  VALUES  OF 
DERIVATIVES 

I 


i  CALL  CNTRL  I 


COMPUTE  VARIABLES 
AND  DERIVATIVES  AT 
TWO  HALF  STEPS 
I 

COMPUTE  VARIABLES 
AT  END  OF  INTERVAL  | 


TEST  FOR  FIXED  OR 
VARIABLE  INTERVAL 

I 

IF  INTERVAL  IS  VARIABLE 
COMPUTE  ERROR.  IF  TOO 
LARGE,  DECREASE  INTERVAL, 
REPEAT  STEP.  IF  TOO  SMALL, 
INCREASE  INTERVAL  AND 
ACCEPT.  _ 


H 


CALL 


CNTRL  (HTHY) 


IF: 

NTRY=1, 

NTRY=2, 

NTKY=3, 

NTRY-U, 

Compute  Next  Step 

Exit  to  Main  Program 
Repeat  Step 

Restart  Integration 

Subroutine  DERIV  Computes 
Derivatives 


Subroutine  CNTRL  Executes  Variable 
Changes,  Stores  Current  Values, 
Executes  Exit  Tests 


This  Loop  Calls  DERIV  0  Times 


NTRY  is  Re -set  in  CNTRL 
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DISCUSSION  OF  RESULTS 


Thia  report  contains  two  methods  for  calculating  the  velocity  potential 
along  sonic  ray  lines  emanating  from  any  point  in  a  non-uniform  flow  field, 
i.e.,  one  that  varies  from  locally  subsonic  to  supersonic  speeds.  Both 
methods  apply  to  pulses  emitted  by  sources  or  doublets.  It  has  been  demon¬ 
strated  that  both  methods  yield  nearly  identical  ray  paths  and  times  of 
transmission.  Those  presented  were  obtained  using  the  second  method. 

Figures  6  through  13  show  ray  paths  of  acoustic  signals  emanating  from 
various  points  in  a  non-uniform  transonic  flow  field.  The  reader  may  want 
to  try  his  hand  at  tracing  one  of  the  ray  paths  in  a  region  of  interest 
such  as  near  a  leading  edge.  If  so,  it  should  be  helpful  to  recall  the 
discussion  starting  with  Equation  (7),  through  the  difference  equations 
of  the  path,  Equation  (8),  and  to  the  end  of  that  section.  An  analysis 
of  the  differential  equation  of  the  path,  Equation  (24)  should  also  be 
helpful.  These  show,  for  instance,  that  where  the  Mach  number  iB  constant 
the  curvature  of  the  ray  path  is  zero;  for  a  given  Mach  number  and  slope 
of  ray  path  the  curvature  is  proportional  to  rate  of  change  of  Mach  number 
along  the  path.  Figures  6,  7>  9,  and  10  conclusively  show  that  when  the 
variation  in  Mach  number  is  parabolic  in  the  chordwise  and  spanwise  direc¬ 
tions  focal  polntsexlst,  both  in  subsonic  and  supersonic  portions  of  the 
flow.  None  of  the  present  theories  accounts  for  the  corresponding  multiple 
crossings  of  the  acoustic  wave  front.  Figures  9  and  12  show  acoustic 
signals  traveling  from  regions  of  supersonic  flow  to  regions  of  subsonic 
flow.  This  can  occur,  of  course,  only  when  the  sonic  line  is  swept  down¬ 
stream.  Figures  9  and  12  also  show  rays  that  have  been  trapped  on  the  Mach 
wave,  travel  outward  to  the  sonic  line  where  the  spanwise  slope  of  the  ray 
path  becomes  zero,  and  are  cancelled  there.  A  study  of  the  ray  paths  that 
cross  the  leading  edge  shows  that  in  practical  applications  it  is  correct 
to  assume  they  do  not  return. 

These  results  permit  the  formulation  of  a  numerical  procedure.  A  box 
method  is  outlined  in  Appendix  III.  It  establishes  velocity  potentials  at 
all  box  centers  on  cm  aerodynamic  surface  and  the  corresponding  generalized 
forces. 
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ACOUSTIC  PATHS 


Figure  c.  Ray  Paths  for  a  Source  or  Doublet  at  (0.l8c,  0.0) 
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ACOUSTIC  PATHS 
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ACOUSTIC  PATHS 


Y 


Figure  o.  Hay  Paths  for  a  Source  or  Doublet  at  (0.6c,  0) 
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ALCUST1L  PATHS 
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ACOlSTIC  PATHS 
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ACOUSTIC  PATHS 
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CONCLUSIONS  AND  RBCOIMENDATION8 


Two  methods  hare  been  outlined  in  detail,  and  one  of  them  has  been 
completely  mechanised  for  calculating  the  velocity  potentials  along 
acoustic  ray  paths  emanating  from  any  point  In  a  non-uniform  transonic 
flow  field  over  a  lifting  surface.  The  one  mechanised  gives  the  ray 
path  and  velocity  potential  for  the  minimum  time  of  travel  from  the 
source  point  to  the  field  point. 

To  calculate  pressures  over  the  planform  and  generalised  forces,  It 
will  be  necessary  to  develop  a  procedure  for  calculating  the  velocity 
potential  at  an  arbitrary  point  due  to  a  sheet  of  sources,  covering  the 
vlng  surface,  and  the  flow  field  In  the  plane  of  the  vlng  out  to  a  dis¬ 
tance  of  several  vlng  spans  In  the  y-dlrectlon,  or  due  to  a  sheet  of 
doublets  covering  the  vlng  surface.  The  latter  is  recomended  for  economy 
reasons. 

The  computer  program  in  this  report  may  be  used  to  refine  the  doublet 
box  method  of  Rodemlch  (3)  in  such  a  way  as  to  Include  the  (possibly  very 
important)  influence  of  vlng  thickness  distribution  on  transonic  airloads. 

A  doublet  box  method  similar  to  the  one  Rodemlch  developed  (Reference  3) 

Is  recooawnded.  The  procedure  Is  heurlstically  described  In  Appendix  III. 
Por  each  of  a  selected  set  of  point  a  In  a  sending  box,  the  distribution  of 
velocity  potentials  along  ray  lines  throughout  the  cone  of  Influence  can 
be  determined.  An  Interpolation  scheme  will  yield  from  these  the  velocity 
potentials  at  box  centers  and  a  numerical  integration  procedure  will  yield 
a  velocity  potential  influence  coefficient  for  each  of  the  box  centers. 

It  vi  11  be  necessary  to  solve  a  set  of  simultaneous  equations  to  establish 
the  strengths  of  doublets  required  to  satisfy  the  tangential  flov  condition 
In  the  subsonic  flov  region.  The  order  of  the  set  will  be  equal  to  the 
number  of  box  centers  In  the  subsonic  region  on  the  vlng.  In  the  super¬ 
sonic  region  the  doublet  strengths  can  be  established  sequentially.  The 
use  of  doublets  to  solve  unsteady  supersonic  flov  problems  has  been  out¬ 
lined  by  Ashley  in  Reference  7. 

It  Is  recoHsended  that  this  method  be  fully  developed  for  the  pur¬ 
pose  of  calculating  generalised  forces  on  vlngs  In  harmonic  motion  at 
transonic  speeds.  A  computer  program  that  would  predict,  vlth  reasonable 
accuracy,  the  flutter  characteri sties  and  unsteady  aerodynamic  loads  on  a 
vlng  responding  to  externally  applied  forces,  such  as  gusts,  would  fill  an 
important  gap  In  available  technology. 
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a  n 


APPENDIX  1.  Prog  ran  Listing* 


sisrrc  main  soo  snicooos 

C  F OR  TRAN  PROGRAM  TO  COMPUTE  (AMO  PLOT)  THE  PATHS  OF  ACOUSTIC  SIC  -  SNICOUIO 
C  NALS  (ANO  TRANSMISSION  TIMES)  ON  AN  AIRFOIL  IN  A  SONIC  FLOW  F  IELD , SN J C001 3 
C  ACCOUNTING  FOR  VARIATION  IN  LOCAL  MACH  NUMBER.  SNIC0020 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CM  =  COEFFICIENTS  OF  MACH  t GUATION.  (SEE  SUBROUTINE  FHACH  )  SNIC002S 
P LX  ANO  PLY  ARE  CONSTANTS  DESCRIBING  THE  PLANFORM  GEOMETRY.  SNIC0030 
THE  PROGRAM  ALLOWS  FOR  EITHER  X  OR  Y  TO  BE  THE  INDEPENDENT  VARIA-  SNICD033 


8LE,  DEPENDING  ON  THE  CURRENT  VALUE  OF  X-P RIME,  WHICH  SETS  IVAR. 


IF  IVAR  =  1, 

YY  =  CURRENT  VALUE  OF  X 
DYYs  CURRENT  VALUE  OF  OX 
XX (1)=  CURRENT  VALUE  OF  Y-PRIMC 
XX 12)  s  CURRENT  VALUE  OF  Y 
XX (3)  =  CURRENT  VALUE  OF  TIME 
XX (4)  =  CURRENT  VALUE  OF  R-OAR 
DXX ( 1)  =  Y-OOUBLE  PRIME 
DXX(2)=  CURR.  VALUE  OF  T-PRIME 
DXX (3) =  CURR.  VALUE  OF  OT/DX 
DXX (4) =  CURRENT  VALUE  OF  DR/OX 

IVAR  IS  ORIGINALLY  SET  IN  MAIN 
PASS  THROUGH  SUBROUTINE  CNTRL. 


IF  IVAR  =  2, 

YY  =  CURRENT  VALUE  OF  Y 
DYY=  CURRENT  VALUE  OF  DY 


SNICD04D 

SNIC0045 

SNICQ090 

SNIC0093 


XX  (1) =  CURRENT  VALUE  OF  X-PRIME  SNIC0060 


XX  (2)  ;  CURRENT  VALUE  OF  X 
XX  (3)  =  CURRENT  VALUE  OF  TIME 
XX  (4)  =  CURRENT  VALUE  OF  R-BAR 
DXX  (1)=  X -DOUBLE  PRIME 
DXX (2) =  CURR.  VALUE  OF  X-PRIME 
DXX  (3) =  CURR.  VALUE  OF  DT/DY 
DXX  (4) 2  CURRENT  VALUE  OF  DR/DY 

PROGRAM,  AND  THEN  RESET  ON  EACH 


WORK  =  WORKING  AREA  FOR  SUDROUTINE  RKS3  . 

IFVD  2  FALSE  ANO  ICKP2  TRUE  FOR  VARIABLE  INTERVAL. 

IFVO  2  true  FOR  FIXED  INTERVAL. 

SX  2  VECTOR  CONTAINING  COMPUTED  X-  VALUES. 

SXP  2  VECTOR  CONTAINING  COMPUTED  X-PRIME  VALUES. 

SY  CONTAINS  COMPUTED  Y  VALUES 
SYP  CONTAINS  CO^:r-UT€0  R-SAR  VALUES 
TIM  CONTAINS  TRANSMISSION  TIMES. 

FM  2  CURRENT  MACH  NUMBER 

ISORS  2  -i  CCFINCS  A  SUPERSONIC  SOURCE,  RECEDING  PATH. 
ISROS  r  0  DEFINES  A  SUPCRSONIC  SOURCE,  ADVANCING  PATH. 
ISORS  =  l  CCFINCS  A  SUOSONIC  SOURCE. 

ior  2  i  for  right  branch,  2  for  left 


)F  X  SNIC0069 

)F  TIME  SNIC0070 
)T  R-BAR  SNICO0TS 
SNIC0000 
X-PRIME  SNICO009 
DT/DY  SNIC0090 
>F  DR/DY  SNICOD93 
SNIC0100 
ON  EACH  SNIC0109 
SNICOl 10 
SNIC0119 
SNIC0120 
SNIC0129 
SNIC013O 
SNIC0I35 
SNIC0140 
SNIC0149 
SNIC0190 
SNIC0199 
SNIC0160 
SNIC0169 
SNIC0170 
SNIC01T3 
SNIC0100 
SNIC0109 
WHEN  NCNTSNIC0190 


ISORS  =  1  DEFINES  A  SUOSONIC  SOURCE.  SNICO100 

IOR  2  1  FOR  RIGHT  BRANCH,  2  FOR  LEFT  SNIC0103 

NCNT  IS  THE  COUNTER  FOR  THE  VECTORS  SX , SY, SXP , SYP , T IM.  WHEN  NCNTSNIC0190 
2  NMAX,  IMTEGfi AT  ION  STOPS,  AND  THE  FLOW  PASSES  TO  NEXT  PATH  SNIC0193 

I  TRAP  21  INDICATES  SIGNAL  IS  TRAPPED  ON  THE  LOCAL  MACH  CONE.  SNIC02OO 

0/  IS  INITIAL  VALUE  OF  INCREMENT.  SNIC0203 

CINF  2  REMOTE  SPEED  CF  SOUNO  IN  ROOT  CHORDS  PER  SECOND.  SNIC0210 

FMINF2  REMOTE  MACH  NUMBER  SNIC0213 

PwTE  -  THE  POTE  MATRIX  CONTAINS  THE  VELOCITY  POTENTIALS  ALONG  A  SNIC0220 
RAY  PATH,  NORMALIZED  ON  BO  .  SNIC0229 

FREQ  = ASSUMED  FREQUENCIES  IN  RADIANS  PER  SECOND.  SNIC0230 

SNIC0239 

EXTERNAL  DCRIV,  CNTRL  SNICO240 

COMMON  SNIC0249 

4/WORK/  WORK (90)  SNICO290 
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o  n 


♦/xyz/  sx  non  ,sxp  non  , sy  non  , syp  non  ,al  (4n  , Tinnon 

♦/XDX/  XX (4)  ,DXX  (4) , YY,DYY,DZ 

♦  /CM/  CM  (6) 

♦  /TABLE/  ATABL  (4)  ,RTA0L (4) 

♦/PL/  PLX (0) ,PLT (0) 

♦  / I CNT/  I  VAR , NCNT , I  SOP  5 ,  I  BR , I  TRAP , NMAX 

♦  /SOURCE/  XO(20)  ,  YO(20) 

♦  /EPS/  El ,£2,FM,  YMAX 
♦/NNN/  N5S , NLC3 , NLL3 
4/ECM/  ECH 

4/C4/  CM2 (7) 

C 

1000  FORMAT (2L12  ) 

1010  FORMAT (6E12. 0) 

1020  FORMAT (6112  ) 

3  READ  (5,1020)  NSORCE,NLA,NPL ,NMAX ,NF 
READ  (5,1000)  FVO.IBW* 

READ  (5,1010)  (XO(I)  ,YO(I)  ,  I  =  1,NS0RCE) 

READ  (5,1010)  (CM  (I) ,1  =  1,6  ) 

READ  (5,1010)  02, El  ,£2, YMAX 

READ (5, 1010) (ATABL(I) , 1=1,4) , (RTABL(I) ,1=1,4  ) 

READ  (5,1010)  (PLX(!)  ,PLT(I)  ,  1=1, NPL  ) 

READ  (5,1010)  CINF,FM1NF,TAU,TSAA 

TAU=MAX.  (T/C),  TSAA  =  TANGENT  OF  SEMI-APEX  ANGLE 
DIMENSION  FREQ(IO) ,  POTE ( 101 , 2, 10  ) 

READ  (5,1010)  (FREQ  ( I)  ,  1  =  1  ,NF  ) 

C 

DIMENSION  XSO(40>  ,TSO(40) 

ECM  =  CINF*SQRT(5.0«FMINF4*2) 

ECM= 1 . O/ECM 

CALL  CONK  (40,NXY,YMAX,YSO,XSO,  IER  ) 

2000  FORMAT (49H0  error  in  subroutine  sonic.  CHECK  MACH  CONSTANTS  ) 

GO  TO  (1,2),  IER 
2  WRITE  (6,2000) 

1  CONTINUE 
NVAR  =4 

c  NVAR  IS  THE  NUMBER  OF  VARIABLES 
CM2 (1 )  =0.3 
CM2  (2)  =0.T 
CM2  (3)  =ATAN(1./T3AA) 

CM2 (4)  =TAU 
CM2  (5)  =1.104>TSAA 
CM2 (6)  =.04 
CM2(T)=FMINF 
C  OEVELOP  LAMOAS 

NL=2* (NLA/2) 

C  THERE  WILL  ACTUALLY  BE  NL  VALUES.  IF  NLA  IS  EVEN,NL=NLA.  BUT  NL= 
C  NLA  -  1  IF  NLA  IS  O0O. 

NL1*NL-1 


SNICL255 

SNIC0260 

SNIC0265 

SNIC02T0 

SNIC02T5 

SNICO20O 

SNICO205 

SNIC0290 

SNIC0295 

SNIC0300 

SNIC0305 

SNIC0310 

SNIC0315 

SNIC0320 

SNIC0325 

SNIC0330 

SNIC0335 

SNIC0340 

SNIC0345 

SNIC0350 

SNIC035J 

SNIC0360 

SNIC0365 

SNJC0370 

SNIC0373 

SNIC0380 

SN ICO305 

SNIC0390 

SNIC0395 

SNIC0400 

SNIC0405 

SNIC0410 

SNIC0415 

SNIC0420 

SNIC0425 

SNIC0430 

SNIC0435 

SNIC0440 

SNIC0445 

SNICL430 

SNIC0455 

SNIC0460 

SNIC0465 

SNIC0470 

SNIC0472 

S  .C047S 

CNICO40O 

5NICO403 

SNIC0490 

SNIC0495 


3t> 


NL2  -NL/2 
XN-  NL2*(NL2*1) 

0G  =  6.2031 0/XN 
AL (1 ) =0. 

00  10  J=3,NL1 ,2 
Xi  =  (J-l)/2 
J1  =  J-l 

AL(J)=AL(J-2) ♦XJ*DC 

10  AL(J1)=-AL(J) 

AL(NL)*  3.14199 

C  SET  UP  GRID  LIMITS 
XU=0. 

XL3J. 

YLs-YNAX 
TR  =  YMAX 

CALL  LIMITI ( YL , YR, XL ,XU) 

DO  600  NS=1 ,NSORC£ 

NSS  =  NS 

CALL  GRAPH(1,42, -NPL ,PL  Y,PLX , 2H  Y,2H  X.15H  ACOUSTIC  PATHS  ) 
XOF=XO(NS) 

YOF  =  YO(NS) 

CALL  GRAPH(0,42,-NXY, YSO.XSO  ) 

NLLS  s  NL 
DO  900  NLC=1,NL 
NLCS  =  NLC 
I  TRAP  s  0 

CALL  FMACH  (XOF,YOF ,FM,FMX,FMY  ) 

TEST1  =FM  -  COS (AL (NLC! ) 

TEST2  =  SIN (AL  (NLC)) 

IF (NLC  .NE.  NL  )  GO  TO  11 
IF  ( YOF  .GT.  0.  )  GO  TO  11 
TEST2  =  -TEST2 

11  IF  (NLC-1 )  14,12,14 

12  IVAR=1 
GO  TO  30 

14  IF(NLC-NL)  10,12,10 
10  IF  (TEST1)  22,20,22 
20  I VAR=2 
GO  TO  30 

22  TEST  =  TEST1/TEST2 
ART  =  ABS (TEST) 

IF(ART-l.O)  20,12,12 

30  CONTINUE 
C  SET  IBR 

F  L  =  AL(NLC> 

IF (NLC-1)  32,31,32 

31  IF ( YOF)  41,41,42 

32  IF(NLC-NL)  36,34,36 
34  IF ( YOF)  42,41,41 

36  IF (FL)  42,42,41 


SNIC0900 

SNIC0305 

SNIC0910 

SNIC0919 

SNIC0920 

SNIC0529 

SNIC0930 

SNIC0535 

SNtC0940 

SNIC 1949 

SNICC390 

SNIC0399 

SNIC0960 

SNIC0969 

SNIC09T0 

SNIC09T9 

5NIC05S0 

SNIC05B9 

SNIC0990 

SNIC0599 

SNIC0600 

SNIC 06 09 

SNIC0610 

SNIC 06 19 

SNIC0620 

SNIC0629 

SNIC0630 

SNIC0635 

SNIC0640 

SNICD649 

5NIC06J0 

SNIC0639 

SNIC0660 

SNIC0669 

SNIC 06 70 

5NIC06~9 

SNICO60O 

SNICO609 

SNIC0690 

SNIC0693 

SNIC0700 

SNIC0709 

SN I COT 10 

SN I COT 19 

SNIC0T20 

SNIC0T29 

SNIC0T30 

SNIC0739 

SNIC0T40 

SNICOT49 
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41  IBR  =  1 

CO  TO  50 

42  IBR  =  2 

50  CONTINUE 

C  SET  ISORS 

csl  =cos<fl) 

RM  = 1 . 0/FM 

IF(FM-I.O)  60,51,51 

51  IF((Fm-1.0)-C2)  52,52,51 

52  CO  TO  (53,54) , IVAR 

53  TPR  =TEST2/TE3T1 

TST=  1 . 0- YPR442*  (FH4#2  -  1.0) 

55  IF (TST-E1 )  500,500,5# 

54  XPR  =  TEST1/TE3T2 

TST=  XPR4*2-(FM4#2-1.0) 

CO  TO  55 

5#  IF (CSL -RM)  60,60,64 
60  I SORS= 1 
CO  TO  70 
64  I SOR03  -1 
CO  TO  70 
60  ISORS  =  O 
C 

70  NCNTs 1 

CO  TO  (00,90) , IVAR 

C  IF  IVAR=1,X  IS  THE  IN0EPCN0ENT  VARIABLE. 

C 

so  rr  =  xof 

IF(TESTl)  01,01,02 
01  orr=-D2 
CO  TO  03 
02  orr  =  02 

03  XX  (1)  =  TEST2/TEST1 
XX  (2)  =ror 
XX (3)  =  0. 

XX (4)  =  0. 

CO  TO  100 

c  IF  IVAR-2,r  IS  THE  INCEPENDENT  VARIABLE. 

90  rr  =  roF 

CO  TO  (91,92)  ,  I8R 

91  orr  =  oz 
CO  TO  93 

92  orr  s  -02 

93  XX  (1)  =  TEST1/TEJT2 
XX  (2) S  XOF 

XX (3)  *  0. 

XX (4 )  3  0. 

100  CALL  RKS3  (OCR I V, CNTRL , XX , DXX , A TASL  , RTAOL , WORK , 
IKP.NTRT, ICRR  ) 

1070  FORMAT (1  HI , 30X » 43H  IVAR  NCNT  ISORS  ICR 


5NIC0750 

SNIC07S9 

SNIC0760 

3NIC0769 

SMJC0770 

SNIC077J 

SNIC0700 

M1C07SS 

5NICO70O 

SNICO705 

SNIC0000 

SNICO00J 

SNICO0I0 

5*  I  CO#  15 

SNIC0020 

SNICO025 

SNICO0J0 

S'(tC0039 

5NIC004O 

SNICO045 

S'lICOMO 

S’IIC0J)9 

S'JICOMO 

£-MC0.a:-5 

S'<  IC0T70 

S.WIC0073 

S'MCO.MO 

SNICOMJ 

s,(ico.»;0 

i’<JCU»>5 
SHICO9J0 
S*(IC0>05 
s-iico;io 
SWICC9I9 
S.NICO920 
SNIC03J5 
S4ICOJ30 
SMI CO 93 5 
SH1C0940 
SNIC094J 
SNICO950 
SNIC0955 
SNIC0960 
SNIC0965 
SMICO970 
SNIC0975 
SNICO900 

rr,Orr,NVAR, IFvO, IBSNICO905 

SNIC0990 

I  TRAP  NLC9  =  )  SNIC0995 
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1060  FORMAT (1H0.27X ,  617  > 

WRITE  16,1070) 

WRITE  (6,1060)  I  VAR , NCNT, I $OR  S,  18R,  I  TRAP , NLCS 

1060  FORMAT ( 22H  ERROR  IN  RKS3,  IERR  =  F4  ) 

IF(IERR)  103,140,103 
103  WRITE  (6, 1060)  IERR 
60  TO  300 

1050  FORMAT(lH-,42X,4HXO  =  E16.8/  43X,4HYO  =  E16.8/  43X,10HMAC)i  NO. 
116.8//  29X , 3 1 H  ACOUSTIC  RAY  PATH  FOR  LAMBDA  =  E16 . 8///17X , 1HX , 
11HY,14X,7hX-PRIME,11X,7HR*BAR  , 12X , 4HT ’ME//  ) 

1040  Format ( 1  m  7x, 3E18.8) 

140  WRITE  (6,1050)  XO  (NS'  ,  YO(f.S)  ,FM,FL 

WRITE  (6, 1040)  (SX  (I)  , ST (I) , SXP ( I)  ,SYP  (I) , TIM ( I )  , 1  =  1  ,NCNT  ) 

CALL  GRAPH  (0,NLC , -NCNT , ST, SX  ) 

CALL  POT  (NF.FREQ.POTE  ) 

1100  FORMAT (1H1,25X,54H  VELOCITY  POTENTIALS  ALONG  A  RAY  PATH  FOR  A 
ICE  AT  ) 

1110  FORMA  T  ( 1  H  -  , 42X , 4HXO  =  E16.8/43X,  4HYO  =  E15.S/  43X,  8HLAMDDA  = 
1.0  //39X , 3 OH AL TERN AT  I NO  REAL  AND  IMAGINARY  ) 

1120  FORMAT  ( 1 M- ,6X , 7 NOMEG A  =E16 .«//  ) 

1090  FORMAT  (1H  6X.6E16.6) 

DO  300  N=1,NF 

IF  (N  .  NE .  1  )  GO  TO  200 

WRITE  (6,1100) 

WRITE  (6,1110)  XO(NS)  ,  YO(NS)  ,FL 
200  WRITE  (6,1120  )  FREO(N) 

WRITE (6, 1090)  ((  POTE(I,K,N)  -K=l,2  ),  1  =  1, NCNT  ) 

300  CONTINUE 
500  CONTINUE 
600  CONTINUE 
GO  TO  3 
END 


SN I  Cl  000 
SNIC1003 
SNTC1 010 
SN I  Cl  01 5 
SNIC1020 
SNIC1023 
SNIC1030 
SN I  Cl  035 
=  ESNIC1040 
17X , SNIC1 045 
SN I  Cl  050 
SN I Cl  055 
SN  I C 1 06  0 
SNIC1065 
SNIC1070 
SNIC1073 
SN I Cl  080 
SN;C1085 
SOURSN IC 1 090 
5NIC1P95 
E16SNIC1 iuO 
*N ICI 105 
SNIC1110 
SN I Cl  1 13 
SNICU20 
SNIC1 125 
SNIC1130 
SN I  Cl  135 
SNIC1 *40 
SNIC1  i-vS 
SN  ICi  1 50 
SNIC1155 
SN  ICI 160 
SN I  Cl  163 
SN ICI  170 
SNIC1175 


36 


start c  ocr i  soo 

SN I  Cl  1 80 

SUBROUTINE  OCR  I V 

SNIC1185 

COMMON 

SNIC1190 

♦/XOX/  XX (4)  ,DXX  (4) , TY ,D YY ,OZ 

SNIC1195 

♦/CM/  CM (6 ) 

SNIC1200 

♦  /  I CN  f /  IVAR.NCNV, JSORS,  I BR , ITRAP.NMAX 

SNIC12Q5 

♦/CPS/  Cl ,E2,F», YMAX 

SNK1210 

♦/NNN/  NSS.NLCS.NLLS 

SNIC1215 

♦/ECM/  CCH 

SNIC1220 

c 

SNIC1 *25 

CO  TO  ( 1  0, 50) , I VAR 

SN I  Cl  230 

c 

X  IS  THE  INDEPENDENT  VARIABLE 

SNIC1235 

io  call  fmach  (yy,xx <2) ,fm,fmx,fmy  > 

SN I Cl  240 

R=XX(1) 

SNIC1245 

OXX  (2) SR 

SNIC1250 

0  =FM4fM  -1.0 

SN  I C 1 2.'  5 

TSI  S  1.0-R4R4B 

SN I C 1 26  0 

A  =FM4>FM  ♦  5.0 

SNIC1265 

SA  s  SORT (A) 

SN I  Cl  270 

IF  (B)  103,103,101 

SN I  Cl  275 

101  BETA  s  SORT  (B) 

SNIC1280 

IF  ( ITRAP  .EQ.  n  GO  TO  104 

SNIC1285 

IF  ( I SORS  .EO.  1)  CO  TO  103 

SN I  Cl  290 

IF  (TSI  .GT.  El)  CO  TO  103 

SNIC1295 

ITRAP  s  l 

SNIC1300 

GO  TO  1 04 

SNIC1305 

103  IFfTSI  .GE.  0.  )  GO  TO  215 

SNIC1310 

ITRAP  s  2 

5NIC1315 

TSI  =  0. 

SNIC1320 

215  RAD  =  SORT  (TSI) 

SNIC1325 

DXX  (4 )  s  RAO 

SNIC1330 

RA0=  1 . 0/ (A4B) 

SNIC1335 

TMlrFM4 (FM442  ♦  1 1 . 0) /B 

SNIC1340 

TM2=  2. OoFM^ <FMv42*8.0> AR442 

5NIC1345 

TM3=  ( (RADfO)  /D)  ♦  (7  .  D#rM4A2  *5 . 0) 

SNIC1350 

TM4=(FM/A)AR^(6.04RM2  ♦l.O) 

SNIC1355 

GC  TO  105 

SNIC1360 

1 04  RBB  =  1.0/  (8^2) 

SN I  Cl  365 

OXX  (4)  S  0. 

SNIC1370 

105  IF  ( I SORS)  11,15,15 

SNIC1375 

11  GO  TO  (12, 13) , IBR 

SNIC1390 

12  IF  ( I  TRAP)  91,91,7 

SNIC1395 

13  IF  ( I  TRAP)  92,92,9 

SNIC1390 

15  GO  TO ( 1 6 , 17) , IBR 

SNIC1395 

16  IF  ( I  TRAP)  92,92,7 

SN I Cl  400 

17  IF  ( I  TRAP)  91,91,9 

SNIC1405 

19  GO  TO  (92,91) ,  IBR 

SNIC1410 

91  IF  (R)  4,3,3 

SNIC1415 

92  IF  (R)  3,3,4 

SNIC1420 

c 

T  IS  THE  INDEPENDENT  VARIABLE 

SNIC1425 
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50  CALL  F MACH  (XX (2)  ,  Y  Y , FH , FHX  ,  FRY  ) 
ft  =  XX  (1) 

Oxx  (2)  :  R 
B  =  FM4FM-1.0 
TSI  =  R4R-B 
c 

A  =  i.O*FM*FM 
SA  =  SORT (A) 

IF  (B  .IT.  0.  )  60  TO  100 

106  BETA  =  SORT (B) 

irtiTRAp  .ca.  i)  60  to  109 
irdSORS  .E8.  1)  60  TO  100 
•(TSI  .6T .  El)  CO  TO  100 
I  TRAP  =  1 
CO  TO  109 

100  IF(TSI  .CE.  0.)  CO  TO  107 
T9I=0. 

(TRAP  s  2 

107  RAO  =  SORT (TSI) 

OXX (4)  =  RAO 
RAB  s  1.0/ (APB) 

TM1=  (FM/B)4  (FM4>*2M1 .0) 0R**3 
TM2-  2.O*FM4(FM**2*0.O)*R 
TM3= (RA04»3/B)P(7.04FMP*2a5.0) 

TM4  =  (FM/A)  *  (R4'42+6 . 0) 

CO  TO  110 

109  OXX  (4)  s  0. 

110  IF(ISORS)  32,60,60 
52  CO  TO  (54,56) , IBR 
34  IF ( I  TRAP  >  1,1,5 
56  IF ( I  TRAP  )  2,2,0 
60  CO  TO  (62,64) , IBR 
62  IF  ( I  TRAP  )  2,2,5 
64  IF ( I  TRAP  )  1,1,0 
68  CO  TO  (2,1) , IBR 

C  FORMULAS  FOR  THE  SECOND  DERI  VS  FO'.lOW 
C 

1  IF  ( AOS (B)  .IE.  l.E-03)  CO  TO  220 
DXX(1)=RAB*(-TM1  ♦  TM2  -  TM3) *FM*  ♦  TM4  *FMX 
DXX  (3)  =  (SAPECM/B) *(FM*XX (1) aRAO) 

CO  TO  100 

2  IF  (ABS (B)  .CT.  l.E-03  )  CO  TO  209 

220  0XX(1)  =  (.3/A> *(2.*R**3*R*9./R)*FMY  ♦  (FM/A) * (R**2 »6 . )  4fMX 
DXX  (3):  (1 .224  75*ECM)*(RM1./R)) 

CO  TO  100 

209  OXX(l):  RA04(-TH1*TM2aTM3)4FMY*TM4  *  FMX 
DXX  (3)  =  (SA4ECM/P) #(FM8XX (1) -RAO) 

CO  TO  tOO 

3  IF  (NLCS  .E«.  NLLS)  CO  TO  4 

OXX (1 )  =  RA84 (TM1 -TM2*TM3) 0FMY  -TM4P  fMX 


SNIC 

1430 

SNIC 

1435 

SNIC 

1440 

SNIC 

1443 

SNIC 

1430 

SNIC 

1455 

SNIC 

1460 

SNIC 

1465 

SNIC 

1470 

SNIC 

1473 

SNIC 

1400 

SNIC 

1485 

SNIC 

1490 

SNIC 

1495 

SNIC 

1300 

SNIC 

1503 

SNIC 

1310 

SNIC 

1513 

SNIC 

1520 

SNIC 

1523 

SNIC 

1530 

SNIC 

1535 

SNIC 

1340 

SNIC 

1545 

SNIC 

1550 

SNIC 

1555 

SNIC 

1560 

SNIC 

1565 

SNIC 

1570 

SNIC 

1375 

SNIC 

1500 

SNIC 

1585 

SNIC 

1590 

SNIC 

1595 

SNIC 

1600 

SNIC 

It  J5 

SNIC 

1610 

SNIC 

1615 

SNIC 

1620 

SNIC 

1625 

SNIC 

1630 

SNIC 

1635 

SNIC 

1640 

SNIC 

1645 

SNIC 

1650 

SNIC 

1655 

S  NIC 

1660 

SNIC 

1665 

SNIC 

1670 

SNIC 

1675 

UO 


DXX ( 3 ) -  (SA4CCM/B) ♦  (FM  ♦  RAO) 

00  TO  100 

4  ir(AB3(B)  .CT.  1 .€-03  )  00  TO  203 

204  CXX(l)::-(.5/AH>(9.4R<'*4*R**2*2.)*FMY  -  (R/A)  *  (6  .  *R**2  ♦  1  . )  *FMX 
OXX  (3)  =  (1 .22475*CCM) 4(1  . *R$R) 

CO  TO  100 

205  OXX  ( 1 ) -  RA0*(TM1-TM2-TM3) 4FMY-TM4*  FHX 
DXX (3)  =  (SA4ECM/B) 4(FM  -  RAO) 

CO  TO  100 

5  0XX(1)::FM*(<FMY/BETA)  ♦FMX) 

OXX (3)  =  (5A4ECM/B) *  FH  *  XX  (1) 

CO  TO  100 

6  0XX(l)=FM*((-FMr/8ETA)  4fMX  > 

DXX (3)=  (SA4ECM/D) *FM*XX ( 1 ) 

CO  TO  100 

7  0XX(1)=  -  (FM4R08) * (FMY^BETAAf  MX  ) 

0XX(3)=  (SA4ECM/B) *FM 

CO  TO  100 

8  OXX <1)=FM*RBB*(-FMY4B£TA*THX> 

0XX(3)  =  (5A4ECM/B) 8f  M 

100  IF(0rr  .LT.  0.  )  CO  TO  31 
OXX (3)  =  ABS (OXX (3) ) 

CO  TO  32 

31  DXX (3)  =  -  1.0* (ABS (OXX (3))) 

OXX (4)  =  -1 . 0* (ABS (DXX (4)  )  ) 

32  RETURN 
END 


SNIC1680 
SNIC1605 
SN I C 1690 
SN I  Cl  695 
SN I  Cl  700 
SN I  Cl  705 
SN I  Cl  71 0 
SN I  Cl  71  5 
SNIC1720 
SN I  Cl  725 
SNIC1730 
SNIC1735 
SN I  Cl  74  0 
SN I  Cl  745 
SN I  Cl  750 
SNIC1755 
SNIC1760 
SNIC1765 
SNIC1770 
SNIC1775 
SNIC170O 
SNIC1785 
SNIC1790 
SNIC1795 
SNIC1BOO 
SNIC1805 
SN I Cl  810 


Ul 


jierrc  cont  sdo 

SNI Cl  81 3 

SUBROUTINE  CNTRL (NTRY) 

SNlCiSZO 

COMMON 

SNIC1823 

*/xrz/  sx  non  ,sxp  non  ,  sy  (ion  ,syp  (ion  ,  al(4d  ,tim(iod 

SNIC1830 

4/XCX/  XX (4) ,DXX (4) , YY,OYY,DZ 

SNIC10J5 

*/CM/  CM  (6) 

SNI Cl  940 

*/ ICNT/  I  VAR i NCNT  , ISORS,  IBR, ITRAP ,NMAX 

SNIC1045 

4/EPS/  El ,C2  M , YMAX 

SNICl 850 

4/NNN/  NSS ,NLCS,NLL3 

SNIC1855 

IF (NCNT  .NE.  1  )  CO  TO  ( 

SNICl 960 

NCO  =  1 

SNIC1063 

IF  (NR  .EQ.  I)  CO  TO  t 

SNI C1870 

NR  s  1 

SNI Ci 87 3 

IF(A8S(DXX(1)*DYY  )  .LE.  .25)  OO  TO  6 

SNICl 880 

4  urr  =  .3*orr 

SNI Cl  095 

if(abs(Dxx(1)*dyv  )  . le .  .25)  co  to  r 

SNICl 090 

GO  TO  4 

SNI Cl  093 

r  NTRY  =  4 

SNICl 900 

RETURN 

SNI Cl  905 

6  IF  (A0S  (XX  ( 1 ) )  .IT.  1.0  )  GO  TO  20 

SNIC1910 

1  NTRY  =4 

SNIC1915 

GO  TO  (2,3)  , I VAR 

SNICl 920 

2  IVAR=2 

SNICl 925 

GO  TO  3 

SNI Cl  930 

3  IVAR=1 

SNIC1935 

c 

SWITCH  VARIABLES, SET  NEW  INITIAL  CONDITIONS 

SNIC1940 

5  SAV  =YY 

SNICl 945 

DYY  =  D YY4XX ( 1) 

SNK'1950 

10  YY  =  XX (2) 

SNICl 953 

XX (1)  =1 .0/XX(l) 

SNICl  960 

XX (2)  =SAV 

SNIC1965 

RETURN 

SNIC19T0 

20  GO  TO  (23,35) , IVAR 

SNICl 975 

c 

STORE  CURRENT  VALUES  WHERE  X  IS  INDEPENDENT  VARIABLE. 

SNICl 900 

23  SX (NCNT)  =  YY 

SN IC1 985 

c 

CHANGE  IBR  WHEN  Y-P RIM  PASSES  THROUGH 

ZERO 

SNICl 990 

IF  (ADS (XX  (I ) )  .CT.  1.0  E-02)  GO  TO  13 

SNICl 995 

IF ( (DXX  (1)*DYY*XX  (1) )  .GE.  0.0)  GO  TO 

15 

SNIC2000 

IF (NCO  .£0.  2)  GO  TO  13 

SNIC2005 

NCO  =  2 

SNIC2D10 

XX (1) =-XX  (1) 

SNIC201 5 

NTRY  =  4 

SNIC2020 

GO  TO  (11,12) ,  IBR 

SNIC2023 

11  IBR  -2 

SNIC2030 

CO  TO  19 

SNIC2035 

12  IBR  =  1 

$NIC2r40 

GO  TO  1* 

SNIC2045 

15  IF (NCO  .NE.  2  )  GO  TO  19 

SN I C2050 

IF (ASS (XX  (1 ) )  .LT.  1.0  €-01)  GO  TO  19 

5NIC2055 

NCO  *  1 

SNIC2060 

U2 


19 

ir  mm  .he.  o.o  >  co  to  27 

SNIC2065 

26 

SXP (NCNT) :  UNCEF 

SNI C2070 

CO  TO  26 

SNIC2075 

27 

SXP(NCNT)  =  1.0/XX(1) 

SN  I C2080 

20 

SY  (NCNT)  =  XX (2) 

SNIC2D85 

S YP  (NCNT)  =  XX(4) 

SNIC2D90 

T IM  (NCNT)  =  XX (3) 

SNIC2095 

CO  TO  50 

5NIC2100 

33 

SX  (NCNT) =XX (2) 

SNIC2105 

TIM  (NCNT)  =  XX (3) 

SNIC2110 

SXP  (NCNT) =XX ( 1 ) 

SNIC2115 

SY  (NCNT) = YT 

SNIC2120 

SYP  (NCNT)  =  XX (4) 

SNIC2125 

50 

CONTINUE 

SNIC2130 

C 

NOW  TEST  FOR  EXIT  CONDITIONS 

SNIC2135 

IF  ( I T R AP  .NE.  2)  CO  TO  51 

SNIC2140 

ITRAP  =  0 

SNIC214S 

NCNT  =  NCNT  -  1 

SNI C2 150 

CO  TO  100 

SNIC2155 

51 

IF ( I TRAP)  60,60,52 

SNIC2160 

52 

TEST  =FM-1.0 

SNIC2165 

IF (TEST)  100,100,53 

SNIC2170 

53 

IF  (TEST-E2  )  100, 100, 6J 

SNIC2175 

60 

IF  (SX  (NCNT) )  100.70,70 

SNIC2180 

70 

IF  (SX  (NCNT) -1.0)  80,100,100 

SNIC2185 

60 

AY  =ADS  (SY(NCNT) ) 

SNIC2190 

IF(AY-YMAX)  105,100,100 

SNIC2195 

105 

IF (NCNT-NMAX)  110,100,100 

SNIC2200 

100 

NTRY  =  2 

SNIC2205 

NR  =  0 

5NIC2210 

RETURN 

SNIC??15 

110 

NCNT  =  NCNT  ♦  1 

SN I CZ220 

RETURN 

SNIC2225 

END 

SNIC2230 

43 


iibftc  macm 

c  master  suor.,  m,  MX,  MT 

SUBROUTINE  FMACM (FX ,F T ,FMS , FMXS ,FMY S) 

COMMON 
#/C4  /  CM2 (7) 

EQUIVALENCE  (A.CM2C1)) , <8,CM2(2>) ,  CAL , CM3 (3) ) ,  (TAU,CM2(4)) 
#  CM2 (5) )  ,  (R1 , CM2 (6) )  ,  (FMINF , CM2  (7) ) 

a  y= abs  (Fy) 

AYY  =  ABS (AK*FX) 

IF (AT  .LE.  AYY)  CO  TO  200 
SK  s  l./  (SORT ( 1 . +AK#AK) ) 

T  = (AV-AYY) #SK 

ioo  call  fmaci  (Fx,ayv,fms,fmxs,fmys) 

CALL  FMAC2  (FX , AY , A, B, AL, TAU, 0 1FM,D1MX ,01MY  ) 

CALL  FMAC2  <FX,AYV,A,B, AL , TAU.D2FM ,D2MX ,D2MY) 

C 

FMS  =  FMS  -0.6*FM INF# (D1FM-D2FM) 

FmxS=  FmxS*FMYS#AK-0.6#FMINF#(01MX-D2MX-AK*D2MY> 

FMYS  =  -0.6#FMINF*01MY*CAY/FY) 

IF  (T  .CE.  Rl)  CO  TO  300 
120  CALL  FMACI  (FX ,FY, SM, SMX, SHY  ) 

ARC  -1.97079#T/RI 
SI  =  SIN (ARC) 

5 MO  =  SI#SI 
FMS= (FMS-SM)*SM0  ♦  SM 
FMXS= (FMXS-SMX) #smo  ♦  SMX 
FMYS= (FMYS-SMY)*SMO  ♦smy 
CO  TO  300 

200  CALL  FMACI (FX.FY, FMS, FMXS, FMYS) 

300  CONTINUE 
RETURN 
END 


SNIC2235 
SNIC2245 
SNIC2240 
SNIC2250 
SNIC2255 
(AK,  SNIC2260 
SNIC2269 
SNIC2275 
SN I C22B0 
SNIC2203 
SNIC2290 
SNIC2293 
SNIC2300 
SNIC2305 
SNIC2310 
SNIC2319 
SNIC2320 
SNIC2329 
SNIC2330 
SNIC2333 
SNIC2340 
SNIC2349 
SNIC2390 
SNIC2399 
SNIC2360 
SNIC2369 
SNIC2370 
SNIC2379 
SNIC23B0 
SNIC2389 
SNIC2390 
SNIC2359 
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IIOTTC  MAC2  SDO 

SUDROUT INC  FMAC2 (X, Y , A , 0 , AL , T AU , DELCP , DCXCP , CD YCP > 

C 

C  SUBROUTINE  COMPUTCS  DELTA  CP 

C 

CS  =  COS(AL) 

CS1  =  1  ./  (SQRTd  .  ♦  (  (1  .  -  A)  ♦♦2)*(CS**2)  )) 

CS2-1 . / (SORT ( 1 . ♦ ( (1 . -fl)**2)*(C5**2) )) 

TA  =  SIN  (AL)  /C S 

TA1=(1.-A)*TA 

TA2=(1.-B)*TA 

EP5=TAU/(2„*3.1415927*A*CS) 

EPS1  =  EPSPC3/CS1 

EPS2  =  EPS*A*CS/( (1 . -B) *C32) 

EOS  =  1.0  -  EPS 
E0S1  =  EPS1  ♦  1.0 
E0S2  s  EPS2  ♦  1.0 
S  =  ABS (X/TA) 

S1=(X-A)/TA1 

S2=(X-B)/TA2 

01=AB3(Y-3) 

Q2=ABS(Y*S) 

Q3=ABS(Y-S1) 

04  =  ABS(YM1> 

05  =  AB3 ( Y-S2) 

06  =AB3(Y*32) 

F AC  =2.*C3/TA 
FAC1=2.*C31/TA1 
EAC2=2. PC52/TA2 

DEL  =-FAC*(01**EPS»Q2**EPS'2.*3**EPS) 

DDX=-FAC*(-1./(01**E0S)  ♦  ! . / (02**cUS>  -2. / (S**EDS) )  ♦  EPS  /TA 
ODT  =  -FAC*  ( 1 . / (Q1 **EOS) ♦! . / (02**EDS)  )  *  EPS 

IF  (SI)  10,10,3 
10  DEL CP r  DEL 
DDXCP=  OOX 
DO  TCP =  DDY 


SNIC2400 

SNIC24D5 

SNIC2410 

SNIC2415 

SNIC2420 

SNIC2425 

SNIC2430 

SNIC2435 

SNIC2440 

SNIC2443 

SNIC2450 

SNIC2435 

SNIC246J 

SNIC2465 

SNIC24/0 

SNIC2473 

SNIC2460 

SNIC2405 

SNIC2490 

SNIC2493 

SNIC2500 

SNIC2503 

SN I C25 10 

SNIC2513 

SNIC2520 

SNIC2525 

SNIC2330 

SNIC2535 

SNIC2540 

5NIC2543 

5NIC2550 

5NIC2555 

SNIC256 0 

SNIC2563 

SNIC2370 

SNIC2575 


CO  TO  30  SNIC250O 

5  0EL1=-FAC1*(1./(03**EPS1) ♦l./(04**EPSl)-2./(Sl**EPSl))  SNIC2503 

D0X1=FAC1*(-1 ./(03**EDS1) ♦  1 . / <04 **EDS 1 ) -2. / (S1**EDS1) )  *  EPS1  /TA1SNIC2590 


DDY1=  FAC1*(1 ./(03**EDS1) d  .  / <04**ED31) ) 

IF  (32)  20,20,30 

20  DCLCP s  OEL»  DELI 
DOXCP=  DOX ♦  00X1 
DOYCP=  DDY*  00Y1 
GO  TO  30 

30  0EL2=-FAC2* (1 . / (Q5**EPS2) ♦  ! . / (06 **EP S2) -2. / (S2**EPS2) ) 
00X2=F  AC2* (-1 . / (03**EDS2) ♦! . /(06**EDS2) -2. / (31**E0S2) ) 
00T2  FAC2*(1./(03**C0S2)  ♦! . / (06**E0S2) ) 

DCLCP  a  DEL*  0EL1  ♦  DELS 
DDXCP*  OOX ♦  00X1  ♦  00X2 


EPS1  SNIC2393 

5NIC2600 
SNIC2603 
SNIC2610 
5NIC2613 
SNIC2620 
SNIC2623 
EPS2  /TA2SNIC2630 
EP 52  SNIC2633 

SNIC2640 
SNIC2S43 


DDT CP-  OD  T ♦  DO T 1  ♦  DDT2 

RETURN 

CNO 


SNIC2650 

SNIC2653 

SNIC2660 


1I0FTC  MAC  1  SOD 

SUBROUTINE:  FMAC1  (FX,FY,FMS,FMXS,FMYS) 

c 

C  SUBROUTINE  COMPUTES  MACH  NO,  MX,  MY . 

C  FX  =  X 

C  F MS  =  MACH  NO. 

c  fmys=  Partial  m  w/resp  to  r 

C  EQ.  FOR  MACH  IS 

C  CM(5)*Y**2>CM(6)*Y**4  ) 

COMMON 
*/CM/  CM (6) 

C 

EQUIVALENCE 

1  (  C  ,CM(D)  ,  <  FMO,CM(2> ) 

2  (A3  ,CM  (5) )  ,  (  A4  ,  CM  (f ) ) 

IF (Fx  .EQ.  0.)  CO  TO  3 
ARC!  =  (-C*FY**2)/FX 
ARC  1  =  -  ABS(ARCl) 

IF  (ABS  (ARC  1 ) 

ARC2 

C 

ARC 3  =  Al*  2.  ♦  A2*FX 

ARC4  =  2.*A3*FY  ♦4.4iA4*FY**3 

EX  =  EXP (ARC  1) 

CO  TO  10 
5  FM3  =  FMO 
FMXS  =  0. 

FM Y5  =  0. 

RETURN 

10  FMS  =  FMO  *EX*  ARC2 

FMXS=  EX* ( (-ARC1/FX)*  ARC 2  ♦ARCS) 
PAUL=  -2. *C*F  Y/FX 
FMYS=  EX* (  PAUL*ARC2  ♦  ARC4) 

RETURN 
END 


SNIC2663 

SNIC26T0 

SNIC26M 

SNIC2680 

FY  =  Y  SNIC2603 

FMXS=  PARTIAL  M  W/RFSP  TO  X  SNIC2690 

SNIC2693 
SN IC2700 
SNIC2705 
SNIC2T10 
SNIC2T13 
SNIC27Z0 
SNIC2723 
(  A2  , CM (4 ) ) , SN IC2730 
SNIC2735 
SNIC2T40 
SNIC2743 
SNIC2750 
SNIC2755 
SNIC2760 
SNIC2765 
SNIC2770 
SNIC2773 
SNIC2780 
SNIC2783 
SNIC2790 
SNIC2795 
SNIC28D0 
SNIC20O5 
SNIC201O 
SNIC2815 
SNIC2820 
SNIC2825 
SNIC203O 
SNIC2035 


M  =  CM(2)  *EXP (-CM(1)*Y**2/X) *(CM(3) *X*CM(4)  *X**2  + 


( A 1  , CM (3) ) 


CE.  50 
A1*FX ♦ A2*FX**2 


,)  CO  TO  3 
♦  A3*f  Y**2 


♦  A4*FY**4 


HOrTC  SON  I  SCO 

SUBROUTINE  SONK(NM,NCR,YM,FY,FX,IER  > 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


=  MAX  NO  OF  X,Y  ALLOWED .  MUST  EQUAL  DIMENSION  Cl'  X,Y,  IN  MAIN 
NCR  =  NO  OF  X,r  ACTUALLY  COMPUTED 
VM  =  MAX.  ALLOWABLE  VALUE  OF  Y 


FX  =  X-VALUES 

FY  =  Y-VALUES 

ICR  =  1  IS  NORMAL  RETURN 

ICR  =  2  INDICATES  AN  ERROR 

CM=  MACH  CONSTANTS  IN  THE  EQUATION  M  =  £XP  (  -CM  ( 1 ) ♦ Y$ 

♦CM(4)*X*X»CM(5)*Y*Y*CM(6)*Y**4)  *CM(2). 

THE  SUBROUTINE  COMPUTES  A  SET  OF  X  AND  Y  VALUES  ON 


*2/X)*(CM(3)  *X 
THE  WING  WHERE 


M=  1 


COMMON 
♦/CM/  CM (6) 

OIMENSION  FX (!) ,FY(1) 

ICR  =  1 
C:CM(I) 

FMO=CM(2) 

A 1  =CM(3) 

A2  =CM(4) 

A3  sCM(3) 

A4  =CM (6) 

FIRST  COMPUTE  X  WHCN  Y=0 
ARC  =  A1 **Z  -4 .+AZ* (FMO-I  . ) 

IF (ARC  .CE.  O.O)  CO  TO  2 
1  IER  =  2 
RETURN 

Z  Fx  ( 1 )  =  (. 5/A2) ♦ (-A1 ♦SORT (ARC) ) 
FY(1)  =  0. 

IF  (FX ( 1 )  .LT.  0.0>  CO  TO  1 
IF  (FX (I)  -IT.  1.0)  CO  TO  4 
FX  m  =  (.5/a 2)  ♦( -A I -SORT (ARC) ) 
IF (FX ( 1 )  .LT.  0.0)  CO  TO  1 
IF (Fx ( 1 )  CE.  1.0)  CO  TO  1 
4  NCR  =  2 
10  NCI :  NCR  -  1 

FX  (NCR) =  FX (NCI) ♦.01 
X=  FX  (NCR) 

R  =  C/X 

0  =  X* ( A 1 ♦ A2*X ) 

TO  =  FY(NC1)R« 

TM1  =  A3-R48 
TM2  =  2 . 4A4 -R#A3 
TM3  =  R*A4 

TM4  :  2.#A4^RR(R*8"2.4A3) 

THJ  :  R*(RAA3-4.#A4) 


SN1C2840 

SMl C2S45 

SNIC2850 

SNIC2855 

SNIC2860 

SNIC2865 

SNIC2870 

SNIC2075 

SNIC2880 

SNIC2883 

SNIC2890 

SNIC2895 

SNIC2900 

SNIC2905 

SNIC2910 

SNIC2915 

SNIC2920 

SNIC2925 

SNIC2930 

SNIC2935 

SNIC2940 

SNIC2945 

5NIC2950 

SNIC2953 

SNIC2960 

SNIC2965 

SNIC29 70 

SNIC2975 

SNIC2980 

5NIC2985 

SNIC299 0 

S*!IC2995 

SNIC3000 

SNIC3005 

SNIC3010 

SNIC3015 

SNIC3020 

SNIC3025 

SNIC3030 

SNIC3035 

SNIC3040 

5NIC3045 

SNIC3050 

SNTC305S 

SNI0060 

SNIC30&3 

SNIC3070 

SNIC3075 

SNIC3080 

SNIC3O05 


»M6  -  R4R4A4 
IMA*  =1 

12  €  T  -CxP  (-R4T0) 

FT  =  L'T  (B»A34TO*A4*TO*TO)  ♦FMO-l  . 
rrr  -€7l  (  TP  1  ♦  TM24 T0-TM34T0442) 
rpr T  -FT*  :TM4*TM5PTO*TM6*TOM2) 

HO  =  -FT/FPT 

IF  (  iFTiFPPi)  .GE.  0.0 )  GO  TO  14 
HO  =  . /?*H> 

14  TO  =TO*HO 

IMA*  =  IMA*  ♦; 

1000  FORMAT (52H0  COMPUTATION  FOR  SONIC 
116.6  ) 

IF ( IMA*  .LT.  10  )  60  TO  10 
WRITE  (6,1000)  HO 
60  TO  1 

10  IF (HO  .6T.  .000 1)  60  TO  12 
Fr(NCR)=  SORT  (TO) 

IF (NCR  .GE .  NM  )  CO  TO  20 
IF (FY (NCR)  .GE.  TM)  CO  TO  20 
IF (F* (NCR)  .GE.  1.0)  CO  TO  20 
NCR  =  NCR  ♦! 

CO  TO  10 
20  RETURN 
END 


SNIC3090 
SNIC3095 
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iibptc  po te 

SUBROUTINE  POT  (NPR ,PR  ,P) 

COMMON 

*/xrz/  sx  ( ion  ,sxp  non  ,sr  iioi)  ,stp  non  ,al (4i) , tim (ion 

•/CM/  CM (6) 

*/ICNT/  IVAR.NCNT, ISORS, I BR , I TR AP  , NMAX 
♦  /SOURCE/  XO(20),TO (20) 

•/EPS/  El ,£2,PM, YMAf 
♦/NNN/  NSS ,NLCS,NLLS 
C 

DIMENSION  PR  (10) ,P (101 ,2,10) 

CON=-. 25/3. 1419* 

XS  sxO(NSS) 
rs  =rO(NSS) 

DO  100  Ns  1 , NCNT 
X=SX (N) 

T=ST(N) 

T  S  TIM (N) 

rbar  s  stp(N) 

10  OO  SO  NTs  1 , NT R 
IP (RBAR)  12,14,16 
12  P (N, 1 ,NP) =0. 

P(N,2,NP)sO. 

CO  TO  SO 

14  r  (N,1,NP) sUNDEP 
P (N , 2 , NP ) sONDEF 
CO  TO  30 

16  IP (RBAR  .LC.  l.E-9)  CO  TO  14 
PACT  =  CON/ROAR 
ARC  s  PR(NP)* T 
CO  s  COS (ARC) 

SI  s  SIN (ARC) 

P(N,l,NP)s  COMPACT 
P(M,2,NP)s  -SI4PACT 
SO  CONTINUT 

i  oo  cc  i nuc 

RETURN 

END 
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BY  RKS3,  RUNGE-KUTTA,  F  4,  V13,  SHARE  D2*ATFRKS3 
(OERIV,CNTRL  ,  ?EAL Y , 0 Y , ATABl  ,  RTABL ,  DEL  TAX  ,X ,  XHALF 
,XZERO,  Y,  YHALF,  YZE  RO, D YHALF, DYZERO.DELTAY, REALX 
,DX ,N, IFVO, IPR^iHTRY, IERR) 


1I9RTC  RK  S3  *  RUNGC-KUTTA,  FORTRAN  IV,  VERSION  13,  SHARE  D2*ATFRKS3 
SU3ROUT INC  RK  S3  CCER  I  V  ,  CN TRL , Y , D Y , A T ACL , R TABL , WORK , X , DX ,N , IFVO 
1  ,  IBKP.NTRY, IERR) 

external  deriv,cntrl 

INTEGER  N , NTR Y ,  IERR 

LOGICAL  IFVO.IBKP 

REAL  Y,OY,ATABL,RTABL,X,OX 

DIMENSION  Y (N)  , C Y  (N) , ATABL (N)  , RTABL (N) 

DIMENSION  WORK (1 ) 

C  DIMENSION  WORK 

CALL  RKINT  (OERIV, CNTRL,Y,0Y, ATABL, RTABL, WORK(l) ,WORK(3) ,WORK(3) 

1  .WORK  (7)  , WORK (9)  .WORK  (2* H*9)  .WORK (4*N»9)  .WORK <6*N*9) 

2  ,WORK(7*N»9) ,WORK(0*N»9) ,X,DX,N,  IFVD, IBKP ,NTRY, IERR) 
RETURN 

END 

SIBFTC  RKINT*  CALLED 
SUBROUTINE  RKINT 

1 
2 

external  ceriv.cntrl 

INTEGER  N,NTR Y, IERR 
LOGICAL  IFVD, IBKP 

REAL  REaLY.OY, ATABL, RTABL, OELTAX.D YHALF, DYZERO.DELTAY, REALX,  OX 
DOUBLE  PRECISION  X,XHALF,XZERO» Y, YHALF , YZERO 

DIMENSION  REALY(N) ,DY(N)  ,  ATABL  <N> , RTABL (N) , Y  (N) .YHALF  (N)  ,  YZERO  (N) 
1  ,DYHALF(N) ,0YZERO(N) .OELTAY(N) 

IERR  =  0 
10  OELTAX  =  DX 
X  =  REALX 
DO  20  1=1, N 
20  Y ( I )  =  REAL Y ( I) 

CALL  OERIV 
GO  TO  200 

30  IF  COX  .Efl.  0.)  GO  TO  230 
OELTAX  =  OX 
0X2  =  OX/2. 

0X4  =  OX/4. 

XZERO  s  X 
DO  40  1  =  1, H 
YZERO  (I  I  =  r  CD 
40  OYZERO(I)  =  0  Y  ( I) 

DO  110  J-1,2 

XHALF  =  X 

X  =  X  *0X4 

REALX  =  X 

DO  30  1=1, N 

OELTAYCI)  =  0 Y  (I) *0X4 

YMALF(I)  =  TCI) 

Y  ( I )  =  YCI)  ♦OELTAY(I) 

50  REAL Y (I)  «  T(I) 
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CALL  OCR  I V 
00  60  1=1, N 

OCLTAYU)  =  CCLTAY  U  >  »D Y (I ) 40X2 
YU)  =  YHALF(I)  ♦DY(I)*DX4 
60  RCALY(t)  =  Y(I) 

CALL  OCR  IV 

X  =  xhalf*oxc 

RCALX  =  X 
DO  70  1  =  1  ,N 

OCLTAYU)  =  OCLTAYU)  ♦DY(I)  *0X1 
YU)  =  YHALFU)  ♦DYU)»0X2 
70  RCAL  Y  (I)  =  YU) 

CALL  OCR IV 
DO  80  1  =  1, N 

DCLTAYU)  s  (OCLTAYU)  *DY ( I ) 00X4) /3 
Y(I)  =  YHALFU)  ♦OCLTAYU) 

60  RCALYU)  =  YU) 

CALL  OCR IV 
GO  TO  (90, 110) i J 
90  DO  tOO  1  =  1, N 
100  DYHALF (I)  =  DY(I) 

110  CONTINUE 

IF  UFVO)  CO  TO  200 
CRRHAX  =  0 
DO  120  1  =  1,  N 

CNR  =  ATACI.  ( I )  ♦ABS  (RTAGL  (I)  0REAL Y  (I)  ) 
IF  (ERR  .CO.  0.)  CO  TO  220 


160  OX  =  OX/1.9040992 

IF  (.NOT.  IEKP)  CO  TO  100 
CRRHAX  =  CRRHAX/10. 

IF  (CRRHAX  .CT.  1.)  CO  TO  160 
CO  TO  100 

170  OX  =  OX/1.9040992 
CO  TO  200 
100  X  =  XZCRO 

DO  190  1  =  1,  N 
YU)  =  Y2CROU) 

190  0 Y  ( I)  =  OYZCROU) 

CO  TO  90 
200  NTRY  «  1 

CALL  CNTRL  (NTRY) 

CO  TO  (30,210,100,10) , NTRY 
210  RETURN 
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SR  =  (DYZCRO(I)  M.fDYHALFU)  ♦CY(I))/3.0OX2 
120  CRRHAX  =  AMAX1 (CRRHAX, A0S  (SR- (RCALY (I ) -SNCL ( YZCRO ( I ) ) ) ) /ERR) 
IF  (CRRMAX-1.)  130,170,160 
130  IF  (CRRHAX-. 73)  140,200,170 
140  IF  (CRRHAX-. 073)  190,200,200 
130  OX  =  0X01.9049992 
CO  TO  200 
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220  I  ERR  =  1 
RETURN 

2 30  1ERR  =  -1 
RETURN 
END 
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FORTRAN  FIXED  10  DIGIT  DECIMAL  DATA 

PECK  NO - PROGRAMMER  _ _  DATE _  PAGE _ of  JOB  WO. 

NUMBER  IDENTIFICATION  DESCRIPTION  DO  NOT  KEY  PUNCH 


FORTRAN  FIXED  10  DIGIT  DECIMAL  DATA 

DECK  NO.  _  PROGRAMMER  .  DATE  PAQE  of  JOB  NO. 

NUMBER  IDENTIFICATION  DESCRIPTION  DO  NOT  KEY  PUNCH 
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FORTRAN  FIXED  10  DIGIT  DECIMAL  DATA 

DECK  fcO. - PROGRAMMER  DATE _  PAGE  of  JOB  NO. 

number  identification  description  do  not  key  punch 
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APPBTDH  III.  Application  to  tha  Boundary  Value  Frcblem 


A  procedure  that,  may  be  used  to  natch  the  tangential  flow  condition 
on  a  wing  surface  is,  In  principle,  the  saae  as  that  employed  by 
Roderalch  in  the  box  method  for  ur  'on  sonic  flow  (Reference  3).  The 
velocity  potential  at  a  field  pci..  (x,y,x)  dun  to  a  doublet  sheet  in 
its  zone  of  influence,  is 


f  r  ^(P,T)^o(x-f,y-T|,z)dfd<n 
S*-V 


(39) 


where  A0(e,v)  is  the  velocity  potential  <5i r continuity  through  the  doublet 
sheet  over  the  region  8  +  W  (the  surface  and  its  wake),  and 


e’ia*n 


(40) 


where  R  -  V  (x-?)2  +  [l-I^U.y,*)]  [(y-T))2  ♦  z2] 

and  where  If  represents  the  number  of  tines  the  wave  front  passes  the 
field  point.  In  uniform  subsonic  f..ow  If  equals  one,  in  uniform  super¬ 
sonic  flow  it  equals  two,  and  in  the  limiting  case  of  uniform  sonic 
flow  it  equals  one.  As  discussed  previously,  in  uniform  sonic  flow  the 
stationary  portion  of  the  perturbation  wave  front  is  not  augmented  by 
high  frequency  signals  that  follow  it;  instead,  the  pressure  discon¬ 
tinuity  is  dissipated  by  them. 

When  the  local  flow  in  a  non-uniform  flow  field  is  sonic  the  wave 
front  gradually  becoaws  stationary  and  is  dissipated.  Rays  of  this 
type  are  ahown  in  Figures  9>  12,  end  13  •  In  certain  regions  of  non- 
uniform  flow  a  wave  front  may  pass  field  points  more  than  twice  as  shown 
in  Figures  6,  7,  9,  10,  and  12.  These  regions  may  be  in  the  region  of 

subsonic  flow  or  in  supersonic  flow.  Multiple  crossings  normally  occur 

on  receding  portions  of  the  wave  front.  Ray  lines  on  advancing  portions 
normally  pass  over  the  trailing  edge  before  they  cross.  In  these 
regions  of  multiple  crossings  of  the  wave  front,  care  must  be  taken  to 
establish  an  accurate  value  of  IT,  and  of  each  of  the  corresponding  gn's, 
n  ■  1,  2,  ...  H.  A  computer  program  that  may  be  used  to  do  this  is 
contained  herein.  Figures  11  and  13  show  that  in  some  regions  of  both 
nub  sonic  and  supersonic  flow  even  tne  receding  ray  lines  do  not  cross. 
All  of  Flguz  -s  6  through  13  show  that  once  a  ray  crosses  the  transition 

region  at  the  edge  of  the  planform  it  does  not  return  to  the  wing 

region.  This  characteristic  is  ljqportcit  because  when  a  doublet  solu¬ 
tion  is  employed  a  ray  trace  can  be  ignored  once  It  reaches  an  edge 
that  Id  not  adjacent  to  the  wake. 
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The  next  step  la  the  procedure  is  to  define  a  grid  of  nquare  boxes 
over  the  region  8  ♦  W,  and  assise  that  , T)  Is  constant  over  the 
area  of  etch  box.  For  this  to  be  a  valid  assuqition  as  nany  as  50  boxes 
along  the  root  chord  nay  be  required.  The  up  wash  adjacent  to  the  upper 
surface  nay  be  written 

W(x,y,0+)  -  Lin  UiiLtll 
z-O*  * 


or, 

l  l  t(x1-f,y1-‘n)dfdT1  (4l) 


l.e.,  the  upvash  at  (x^y^)  equals  the  sunnatlon  (over  all  boxes  B^j, 

that  Influence  it),  of  products  of  the  constant  velocity  potential  dis¬ 
continuities  and  their  downvash  Influence  coefficients.  The  latter  are 
represented  by  the  double  Integral  of  the  kernel  t  over  the  areas  of  the 
boxes.  The  Halts  of  Integration  and  of  Equation  (39)  ere  not 
functions  of  z,  so  frost  Equation  (40)  we  get 


t(Xi'',3rj’7') 


Ad.  i|- 

2VOf  1  *z 


T 


e"lu«n 


(42) 


At  this  point  It  Is  theorized  that  for  non-uniform  flow  around  a  nearly 
planar  surface  the  variation  in  signal  transmission  time  with  distance 
normal  to  the  surface  It  approximately  equal  to  the  variation  In  uniform 
f low,  l.e., 

.  2.  "U-0  1 1 

41  41  C(lf-l) 


or,  performing  the  differentiation 


*n 


±z 

CK 


(43) 


where  the  upper  sign  refers  to  the  advancing  portion  of  the  wave  front  and 
the  lover  sign  to  the  receding  portion.  C  Is  the  speed  of  sound.  Making 
use  of  equation  (43)  when  taking  the  derivative  in  equation  (42). 


♦(*!-'*/! -Tl) 


^i  e2c±ia*R  \  -ioogQ 

2tt  cr3  ^ 


(44) 


The  gn,s  are  those  obtained  by  tracing  ray  paths  through  the  non-uniform 
flow  field. 


One  way  in  which  Equation  (44)  nay  be  evaluated  and  integrated  la  as 
follows:  Say  for  nine  values  of  ( F , Tj )  on  each  sending  box,  the  values 
of  the  kernel  at  the  center  of  the  receiving  box  (x^>/j)  are  evaluated. 

Since  the  ray  paths  are  not  known  in  advance,  each  of  these  values  must 
be  interpolated  from  valuer  in  its  neighborhood.  It  is  then  necessary 
to  evaluate  the  integral  in  Equation  (4l)  given  the  values  of  the  inte¬ 
grand  at  nine  points  in  the  region  of  integration. 

The  unknowns  in  Equation  (4l)  are  the  's.  When  the  center  of 

a  receiving  box  (x1>Jrj)  lles  in  the  subsonic  flow  region  it  lies  in  the 

zone  of  influence  of  every  other  point  in  the  subsonic  region  and  may  lie 
in  the  zone  cf  influence  of  a  small  portion  of  the  supersonic  region 
(Figure  9)-  'HI  velocity  potentials  in  zones  of  mutual  Influence  must 
be  determined  simultaneously.  Once  velocity  potentials  have  been  es¬ 
tablished  that  meet  the  tangential  flow  conditions  on  the  surface  and  the 
zeio  pressure  difference  condition  on  the  wake  they  may  be  fitted  with 
analytical  expressions  that  have  the  proper  edge  behavior.  Using  these 
expressions,  local  oscillatory  pressures  and  generalized  forces  may  be 
obtained  in  the  way  outlined  in  Reference  3 . 
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